rm(list = ls())

wd <- ".../Replication/"
setwd(wd)

# Load/install packages --
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  foreign, 
  ggplot2, 
  estimatr,
  texreg,
  xtable,
  fastDummies,
  sandwich,
  dplyr,
  janitor, 
  gridExtra,
  gsheet,
  zoo,
  interflex,
  lubridate,
  tidyverse,
  stringi,
  readxl,
  ri2,
  modelsummary,
  ggpubr
)

options(scipen=999)

# Note: given the random component inherent to the randomization tests performed
# by ri(), you may observe small differences in the p-values of the randomization
# tests performed using that function.

# Load data ----
d <- read.csv("Data/baseline database.csv")
data <- subset(d, d$group!="Liberado por eleccion")

data$imprisoned_harvest <- (data$dec+data$jan+data$feb)>0
data$imprisoned_noharvest <- (data$mar+data$apr+data$may+data$jun+data$jul+data$aug+data$sep+data$oct+data$nov)>0

declare <- declare_ra(N = nrow(data), m=sum(data$lottery_winner)) 

# Harvest months
harvest<-summary(conduct_ri(formula = imprisoned_harvest ~ lottery_winner,
                            declaration = declare, sharp_hypothesis = 0, assignment = "lottery_winner",
                            data = data, sims = 10000))

with(data, difference_in_means(imprisoned_harvest ~ lottery_winner))

# Non-harvest months
non_harvest<-summary(conduct_ri(formula = imprisoned_noharvest ~ lottery_winner,
                                declaration = declare, sharp_hypothesis = 0, assignment = "lottery_winner",
                                data = data, sims = 10000))

with(data, difference_in_means(imprisoned_noharvest ~ lottery_winner))


Table_4<-data.frame(matrix(nrow = 4, ncol = 3))
colnames(Table_4)<-c('','Harvest Months','Nonharvest Months')

Table_4[1,1]<-'Difference in Means'
Table_4[2,1]<-''
Table_4[3,1]<-'ri p-value'
Table_4[4,1]<-'Num.Obs.'

Table_4[1,2]<-harvest[,'estimate']
Table_4[2,2]<-''
Table_4[3,2]<-harvest[,'two_tailed_p_value']
Table_4[4,2]<-nrow(data)

Table_4[1,3]<-non_harvest[,'estimate']
Table_4[2,3]<-''
Table_4[3,3]<-non_harvest[,'two_tailed_p_value']
Table_4[4,3]<-nrow(data)

View(Table_4)
print(xtable(Table_4), file = "Output/Table 4.tex")

